clear
graph drop _all
set scheme s1mono

cap cd "INSERT FOLDER PATH HERE"

use "BFS Data - Study 1.dta"
gen Study = 1

append using "BFS Data - Study 2.dta"
replace Study = 2 if Study == .

drop if agentid == .

*Assign Study 2 to Wave 10 (was conducted in a single wave)
replace Wave = 10 if Study == 2

*Label main treatment variable
label define AltTreatmentCode 1 "Control" 2 "Restricted" 3 "Agency" 4 ///
	"WP - Early" 5 "WP - Delayed" 6 "Restricted 2" 7 "Restricted WP"
	
*Define panel
xtset agentid Week

cap log close
log using "BFS Results.log", replace

*** Section 1 ***

*Percentage overall effect of full behavioral food subsidy
xtreg TotalTFV i.AltTreatmentCode i.Study Baseline_TotalTFV i.Wave, re ///
	cl(agentid)
nlcom (effect: _b[4.AltTreatmentCode]/_b[_cons])

*Percentage effect of full behavioral subsidy relative to subsidy alone
nlcom (effect: _b[4.AltTreatmentCode]/_b[2.AltTreatmentCode] - 1)

*Percentage overall effect of regular subsidies: see Table 3 code 
nlcom (effect: _b[2.AltTreatmentCode]/_b[_cons])

*Fraction reporting a desire to increase FV consumption 
	*A number of people completed most of the baseline survey without making 
	*it to the end and being assigned a treatment.  We drop them here first.
drop if BaselineDate == .
drop if BaselineTreatment == "" & AltTreatmentCode == .
egen id_tag = tag(agentid)
tab MoreTFV if id_tag == 1

*Fraction of sample choosing FV Subsidy 
	*First need to create indicator variable for agency
gen choice = (AltTreatmentCode > 2 & AltTreatmentCode < 6) if ///
	AltTreatmentCode != .
sum TFV if choice == 1

*Percentage effect of subsidy with agency relative to subsidy without
	*First need to create indicator variables for subsidy and waiting period
	*in order to estimate regression this comes from (Table 3)
gen Subsidy = (AltTreatmentCode > 1) if AltTreatmentCode != . 
gen wp = (AltTreatmentCode > 3 & AltTreatmentCode < 6 | ///
	AltTreatmentCode == 7) if AltTreatmentCode != .
xtreg TotalTFV Subsidy choice wp Baseline_TotalTFV i.Study i.Wave, re ///
	cl(agentid)
nlcom (effect: _b[choice]/_b[Subsidy])

*Percentage effect of subsidy with waiting period relative to subsidy without
nlcom (effect: _b[wp]/_b[Subsidy])
	
*Percentage effect of subsidy with early choice waiting period relative to 
*subsidy with delayed choice waiting period
	*First need to create indicator variables for delayed and early choices
gen early = (AltTreatmentCode == 4)
gen delayed = (AltTreatmentCode == 5)
xtreg TotalTFV delayed early Baseline_TotalTFV i.Wave if ///
	AltTreatmentCode == 2 | AltTreatmentCode == 4 | ///
	AltTreatmentCode == 5, re cl(agentid)
nlcom (effect: _b[early]/_b[delayed])
	
*** Section 2.4 ***

*Subjects completing Part 1 baseline
sum agentid if Study == 1 & id_tag == 1

*Subjects completing Part 2 baseline
sum agentid if Study == 2 & id_tag == 1

*Percentage reporting any food insecurity at baseline
	*First need to make food security numerical
gen insecure_any = (FoodSecure != "Almost never") if FoodSecure != ""
sum insecure_any if id_tag == 1

*Fraction reporting a desire to increase FV consumption
tab MoreTFV if id_tag == 1

*Fraction reporting SNAP participation	
tab SNAP if id_tag == 1

*Fraction male
tab Male if id_tag == 1

*Fraction outside medium/large urban areas
	*First need to create relevant indicator
gen urban = (CitySize > 2) if CitySize != .
tab urban if id_tag == 1

*** Section 4 ***

*Treatment FV spending means
sum TotalTFV if AltTreatmentCode == 1 
sum TotalTFV if AltTreatmentCode == 2
sum TotalTFV if AltTreatmentCode == 3 
sum TotalTFV if AltTreatmentCode == 5  
sum TotalTFV if AltTreatmentCode == 4
sum TotalTFV if AltTreatmentCode == 6  
sum TotalTFV if AltTreatmentCode == 7 

*** Section 4.1 ***

*Table 2, Panel A
	*First need to create a variable for number of Weeks a shopper completes
egen MaxWeek = max(Week), by(agentid)
sum TotalTFV if AltTreatmentCode == 1
scalar c_sd = `r(sd)'
xtreg TotalTFV Subsidy i.Study , re cl(agentid)
xtreg TotalTFV Subsidy i.Study Baseline_TotalTFV, re cl(agentid)
xtreg TotalTFV Subsidy i.Study Baseline_TotalTFV i.Wave, re cl(agentid)
est sto T2PAC3
reg TotalTFV Subsidy i.Study Baseline_TotalTFV i.Wave if Week==1, cl(agentid)
xtreg TotalTFV Subsidy i.Study Baseline_TotalTFV i.Wave if MaxWeek==4, re ///
	cl(agentid)

*Percentage effect size using column (3)
est restore T2PAC3
disp _b[Subsidy]/_b[_cons]

*Percentage effect size using column (3)
disp _b[Subsidy]/c_sd

*Table 2, Panel B
	*First need to create fractional FV measures
gen TotalFood = ReceiptTotal - NonFood
replace TotalFood = TotalTFV + TotalBG  if NonFood == . | TotalFood < 0 | ///
	TotalFood < (TotalTFV + TotalBG)
gen TFV_Ratio = TotalTFV/TotalFood
gen TotalFood_Baseline = Baseline_ReceiptTotal-Baseline_NonFood
replace TotalFood_Baseline = Baseline_TotalTFV + Baseline_TotalBG if ///
	Baseline_NonFood == . | TotalFood_Baseline < 0 | TotalFood_Baseline < ///
	(Baseline_TotalTFV + Baseline_TotalBG)
gen TFV_Ratio_Baseline = Baseline_TotalTFV/TotalFood_Baseline
sum TFV_Ratio if AltTreatmentCode == 1, d 
xtreg TFV_Ratio Subsidy i.Study, re cl(agentid)
xtreg TFV_Ratio Subsidy i.Study TFV_Ratio_Baseline, re cl(agentid)
xtreg TFV_Ratio Subsidy i.Study TFV_Ratio_Baseline i.Wave, re cl(agentid)
reg TFV_Ratio Subsidy i.Study TFV_Ratio_Baseline i.Wave if Week==1, cl(agentid)
xtreg TFV_Ratio Subsidy i.Study TFV_Ratio_Baseline i.Wave if MaxWeek==4, re ///
	cl(agentid)
	
*Test of FV distribution by subsidy 
ksmirnov TotalTFV, by(Subsidy) exact

*** Section 4.2 ***

*Figure 1, Panel A
	*First, create a new treatment code with a more intuitive order for the 
	*figure.  Also need to create treatment means, SDs, and SEs for the figure
gen NewTreatment = AltTreatmentCode
recode NewTreatment (4=5) (5=4)
foreach var in TFV BG {
egen Mean_Total`var' = mean(Total`var'), by(NewTreatment)
egen SD_Total`var' = sd(Total`var'), by(NewTreatment)
egen N_Total`var' = count(1), by(NewTreatment)
gen TopBar_`var' = Mean_Total`var'+1.96*SD_Total`var'/sqrt(N_Total`var')
gen BottomBar_`var' = Mean_Total`var'-1.96*SD_Total`var'/sqrt(N_Total`var')
}
twoway bar Mean_TotalTFV NewTreatment if NewTreatment < 6, barwidth(0.85) ///
	lcolor(gs6) fcolor(gs6) || ///
	rcap TopBar_TFV BottomBar_TFV NewTreatment if NewTreatment < 6, ///
		lcolor(black) ///
	xtitle("") xlabel(1 "Control" 2 "Restricted" 3 "Agency" ///
		4 `" "Delayed" "Choice" "' 5 `" "Early" "Choice" "', ///
		labsize(medlarge) nogrid valuelabels)  ylabel(0 "$0"  4 "$4"  8 "$8" ///
		12 "$12", angle(horizontal) labsize(medlarge)) ytick(0(2)14) ///
	ytitle("", size(medlarge)) xsize(8) legend(off) name(Means_TFV, replace) ///
	subtitle(Panel A: Mean Spending on Fruits and Vegetables, ///
		size(medlarge)) ///
	nodraw
	
*Figure 1, Panel B
twoway bar Mean_TotalBG NewTreatment if NewTreatment < 6, barwidth(0.85) ///
	lcolor(gs6) fcolor(gs6) || ///
	rcap TopBar_BG BottomBar_BG NewTreatment if NewTreatment < 6, ///
		lcolor(black) ///
	xtitle("") xlabel(1 "Control" 2 "Restricted" 3 "Agency" ///
	4 `" "Delayed" "Choice" "' 5 `" "Early" "Choice" "', ///
	labsize(medlarge) nogrid)  ylabel(0 "$0"  4 "$4"  8 "$8" 12 "$12", ///
	angle(horizontal) labsize(medlarge) ) ytitle("", size(medlarge)) ///
	ytick(0(2)14) xsize(8) name(Means_BG, replace) ///
	subtitle(Panel B: Mean Spending on Baked Goods, size(medlarge)) ///
	legend(order(1 "Mean" 2 "95% C.I.") position(12) rows(1) ring(0) ///
	bmargin(medium) region(col(white)) size(medlarge)) nodraw

*Figure 1, Panel C
xtreg TotalTFV i.NewTreatment Baseline_TotalTFV i.Wave if Study == 1, re ///
	cl(agentid)
coefplot, keep(2.NewTreatment 3.NewTreatment 4.NewTreatment 5.NewTreatment) ///
	vertical levels(99 95 90) msymbol(d) mcolor(white) msize(medlarge) ///
	ciopts(lwidth(8 ..) lcolor(gs12 gs8 gs4)) yline(0, lcolor(gs2)) ///
	graphregion(color(white)) legend(off) ///
	coeflabels(2.NewTreatment = "Restricted" 3.NewTreatment = "Agency" ///
	4.NewTreatment =`" "Delayed" "Choice" "' ///
	5.NewTreatment = `" "Early" "Choice" "', labsize(medlarge)) ///
	ytitle("", size(medlarge)) ylabel(0 "0%" 4.03 "+100%" 8.06 "+200%", ///
	labsize(medlarge)) ytick(-2.015(2.015)10.075) ///
	subtitle(Panel C: Relative-to-Control FV Treatment Effects, ///
	size(medlarge)) name(coeffs_FV, replace) nodraw

*Figure 1, Panel D
xtreg TotalBG i.NewTreatment Baseline_TotalBG i.Wave if Study == 1, re ///
	cl(agentid)
coefplot, keep(2.NewTreatment 3.NewTreatment 4.NewTreatment 5.NewTreatment) ///
	vertical levels(99 95 90) msymbol(d) mcolor(white) msize(medlarge) ///
	ciopts(lwidth(8 ..) lcolor(gs12 gs8 gs4)) yline(0, lcolor(gs2)) ///
	graphregion(color(white)) legend(order(1 "99% CI" 2 "95% CI" 3 "90% CI") /// 
	position(12) rows(1) ring(0) bmargin(large) region(col(white)) ///
	size(medlarge)) coeflabels(2.NewTreatment = "Restricted" ///
	3.NewTreatment = "Agency" 4.NewTreatment = `" "Delayed" "Choice" "' ///
	5.NewTreatment = `" "Early" "Choice" "', labsize(medlarge)) xsize(8) ///
	ytitle("", size(medlarge)) ylabel(0 "0%" 8 "+100%" 16 "+200%", ///
	labsize(medlarge)) ytick(-4(4)20) ///
	subtitle(Panel D: Relative-to-Control BG Treatment Effects, ///
	size(medlarge)) name(coeffs_BG, replace) nodraw

*Figure 1
graph combine Means_TFV Means_BG coeffs_FV coeffs_BG, cols(2) rows(2) ///
	iscale(0.65) name(Figure1, replace)
graph display Figure1, xsize(8.5) ysize(5.5) 
graph export Figure1.pdf, replace

*Unadjusted mean difference between control and early-choice waiting period
sum TotalTFV if AltTreatmentCode == 1 
scalar c_mean = `r(mean)'
sum TotalTFV if AltTreatmentCode == 4
scalar ecwp_mean = `r(mean)'
disp ecwp_mean - c_mean

*Percentage overall effect of full behavioral food subsidy
xtreg TotalTFV i.AltTreatmentCode i.Study Baseline_TotalTFV i.Wave, re ///
	cl(agentid)
nlcom (effect: _b[4.AltTreatmentCode]/_b[_cons])

*Percentage effect of full behavioral subsidy relative to subsidy alone
nlcom (effect: _b[4.AltTreatmentCode]/_b[2.AltTreatmentCode] - 1)

*Overall FV subsidy selection
sum TFV if choice == 1

*Treatment subsidy selection rates
sum TFV if AltTreatmentCode == 3
sum TFV if AltTreatmentCode == 5
sum TFV if AltTreatmentCode == 4

*Subsidy selection rate tests
xtreg TFV i.AltTreatmentCode if choice == 1, cl(agentid) re
test  3.AltTreatmentCode = 4.AltTreatmentCode = 0

*Table 3
sum TotalTFV if AltTreatmentCode == 1
xtreg TotalTFV Subsidy choice wp i.Study, re cl(agentid)
xtreg TotalTFV Subsidy choice wp Baseline_TotalTFV i.Study, re cl(agentid)
xtreg TotalTFV Subsidy choice wp Baseline_TotalTFV i.Study i.Wave, re ///
	cl(agentid)
est store T3C3
reg TotalTFV Subsidy choice wp Baseline_TotalTFV i.Study i.Wave if ///
	Week==1, cl(agentid)
xtreg TotalTFV Subsidy choice wp Baseline_TotalTFV i.Study i.Wave if ///
	MaxWeek==4, re cl(agentid)
	
*Percentage effect of subsidy with agency relative to subsidy without
est restore T3C3
disp _b[choice]/_b[Subsidy]

*Rejectable negative effect size of agency subisdy relative to restricted
test choice = -0.64
test choice = -0.65

*Percentage effect of subsidy with waiting period relative to subsidy without
est restore T3C3
nlcom (effect: _b[wp]/_b[Subsidy])

*WP Effect without agency from Part 2 only
xtreg TotalTFV wp if Study == 2, re cl(agentid)

*Baseline to endline effect of wp
	*First need to generate difference variable
gen diff_TFV = Endline_TotalTFV - Baseline_TotalTFV
reg diff_TFV wp i.Wave if id_tag == 1 & Subsidy == 1 & Study == 1, r

*Table 4
xtreg TotalTFV delayed early Baseline_TotalTFV i.Wave if ///
	AltTreatmentCode == 2 | AltTreatmentCode == 4 | ///
	AltTreatmentCode == 5, re cl(agentid)
est store T4C1
test delayed = early
reg TotalTFV delayed early Baseline_TotalTFV i.Wave if ///
	(AltTreatmentCode == 2 | AltTreatmentCode == 4 | ///
	AltTreatmentCode == 5) & Week == 1, cl(agentid)
test delayed = early
xtreg TotalTFV delayed early Baseline_TotalTFV i.Wave if ///
	(AltTreatmentCode == 2 | AltTreatmentCode == 4 | ///
	AltTreatmentCode == 5) & MaxWeek == 4, re cl(agentid)
test delayed = early

*Percentage effect of subsidy with delayed waiting period relative to 
*restricted subsidy
est restore T4C1
nlcom (effect: _b[delayed]/_b[_cons])

*Percentage effect of subsidy with early waiting period relative to 
*restricted subsidy
nlcom (effect: _b[early]/_b[_cons])

*** Section 4.3 ***	

*Percentage choosing smaller FV subsidy over larger BG subsidy
	*First, need to create reimbursement amount (real and hypothetical), the
	*the amount gained (or lost) from choosing the FV subsidy, and then an 
	*indicator for whether someone choosing the FV subsidy lost money by doing
	*so. We refer to this (casually, not formally) as "commitment".
	gen TFVReimbursement = 0.3*TotalTFV
replace TFVReimbursement = 10 if TFVReimbursement>10
gen BGReimbursement = 0.3*TotalBG
replace BGReimbursement = 10 if BGReimbursement>10
gen GainFromTFV = TFVReimbursement - BGReimbursement
gen TFVCommitment = (GainFromTFV < 0) if choice == 1 & TFV == 1 & ///
	GainFromTFV != .
tab TFVCommitment

*Money sacrificed by TFV committers
sum GainFromTFV if TFVCommitment == 1

*Money sacrificed by TFV committers as a percentage of mean subsidy paid
	*First need to create mean actual subsidy paid
gen sub_paid = TFVReimbursement if TFV == 1
replace sub_paid = BGReimbursement if TFV == 0
sum GainFromTFV if TFVCommitment == 1
scalar gain_mean = `r(mean)'
sum sub_paid 
scalar sub_mean = `r(mean)'
disp gain_mean/sub_mean

*Percentage choosing smaller BG subsidy over larger FV subsidy
	*First, create BG corresponding BG commitment variable
gen BGCommitment = (GainFromTFV > 0) if choice == 1 & TFV == 0 & ///
	GainFromTFV != .
tab BGCommitment

*Money sacrificed by TFV committers
sum GainFromTFV if BGCommitment == 1

*Money sacrificed by BG committers as a percentage of mean subsidy paid
sum GainFromTFV if BGCommitment == 1
scalar gain_mean = `r(mean)'
disp gain_mean/sub_mean

*FV subsidy shows more commitment than BG subsidy
	*First, calculate the size of gains from subsidy choice (commitments are 
	*negative of this), then an indicator for whether that gain is positive (
	*not committed) 
gen GainFromChoice = TFVReimbursement - BGReimbursement if TFV == 1 & ///
	choice == 1
replace GainFromChoice = BGReimbursement - TFVReimbursement if TFV == 0 & ///
	choice == 1
gen GainedFromChoice = (GainFromChoice >= 0) if GainFromChoice != . & ///
	choice == 1
xtreg GainedFromChoice TFV, re cl(agentid)

*Average gain from FV subsidy selection
sum GainFromChoice if TFV == 1

*Average gain from BG subsidy selection
sum GainFromChoice if TFV == 0

*Difference in gains by choice
xtreg GainFromChoice TFV if choice == 1, re cl(agentid)

*Figure 2
	*First, need to bin data for presentability
gen XBin = .
gen YBin = .
forvalues i= 1/40 {
    replace XBin =`i' if BGReimbursement >= 0.25*(`i'-1) & ///
		BGReimbursement < 0.25*`i'
    replace YBin =`i' if TFVReimbursement >= 0.25*(`i'-1) & ///
		TFVReimbursement < 0.25*`i'
}
replace XBin = 40 if BGReimbursement == 10
replace YBin = 40 if TFVReimbursement == 10
egen GroupSize = count(1), by(XBin YBin)
twoway line YBin YBin if TFV == 1 & choice == 1, sort lwidth(medthick) ///
	lcolor(gs10) lpattern(dash) || ///
	scatter YBin XBin [fw = GroupSize] if TFV == 1 & choice == 1, ///
		msize(vsmall) mcolor(gs4) msymbol(Oh) ///
	scheme(s2mono) graphregion(color(white)) legend(off) ///
	ytitle("Subsidy payment from Fruits and Vegetables") ///
	xtitle("(Counterfactual) subsidy payment from Baked Goods") ///
	xlabel(9 "$2" 17 "$4" 25 "$6" 33 "$8" 41 "$10") ///
	ylabel(9 "$2" 17 "$4" 25 "$6" 33 "$8" 41 "$10", angle(horizontal)) ///
	name(Figure2, replace)
graph display Figure2, xsize(8.5) ysize(5.5) 
graph export Figure2.pdf, replace

log close
window manage close graph
window manage close graph
